            
            Info << "USnGrad size: " << U.boundaryField()[patchI].size() << endl;
            vectorField USnGrad(U.boundaryField()[patchI].size());
            forAll(USnGrad, faceI)
            {
                vector y1 = {mesh.boundary()[patchI].Cf()[faceI][0], 
                             mesh.boundary()[patchI].Cf()[faceI][1], 
                             mesh.boundary()[patchI].Cf()[faceI][2]};
                vector y2 = {mesh.C()[mesh.boundary()[patchI].faceCells()[faceI]][0], 
                             mesh.C()[mesh.boundary()[patchI].faceCells()[faceI]][1], 
                             mesh.C()[mesh.boundary()[patchI].faceCells()[faceI]][2]};
                vector deltay = y2 - y1;
                vector u1 = U.boundaryField()[patchI][faceI];
                vector u2 = U[mesh.boundary()[patchI].faceCells()[faceI]];
                vector deltau = u2 - u1;
                USnGrad[faceI] = {Foam::sqrt(Foam::sqr(deltau[0] / (deltay[1]+ SMALL))),   
                                  Foam::sqrt(Foam::sqr(deltau[1] / (deltay[1]+ SMALL))), 
                                  Foam::sqrt(Foam::sqr(deltau[2] / (deltay[1]+ SMALL)))};
            }

            Info << "First: " << USnGrad[0][0] << endl;
